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AN  ESTIMATION  METHOD  OF  VERTICAL  DIFFUSION  PARAMETERS 
IN  THE  MESOSCALE  RANGE 


* 

Lei  Hsiao-en  and  Jen  Chen-hai 

(Institute  of  Atmospheric  Physics,  Academia  Sinica) 


Baaed  oa  the  hypothesis  of  similarity  for  local  turbulent  characteristics  in 
planetary  boundary  layers  and  the  principle  of  dimensional  analysis,  and  using 
statistical  form  of  the  vertical  vortex  diffusion  coefficient  K  as  well  as  Ekman 
spiral  wind  profile,  at  vertical  dispersion  pattern  in  the  100  km  range  is 
derived  by  numerical  analysis  method.  It  is  applicable  to  both  flat  uniform 
land  and  complex  terrain.  The  results  so  far  obtained  agree  fairly  with 
experimental  data. 


I .  INTRODUCTION 

During  the  past  decade,  the  long  range  pollution  arising  from 
cities  and  super-high  chimneys  has  become  worse  and  worse.  This 
has  forced  air  pollution  meteorologists  to  take  an  interest  in  the 
study  of  laws  governing  transport-diffusion  in  the  mesoscale  range 
(10-100  km) .  However,  owing  to  limitations  on  methods  of  observa¬ 
tion  and  on  the  understanding  of  mesoscale  processes,  both  the  diff 
usion  theories  and  field  observations  and  tests  are  still  in  an 
embryo  stage,  and  are  far  from  being  able  to  meet  the  demands  of 
the  various  environmental  problems. 

In  this  paper,  we  attempt  to  use  the  hypothesis  of  similarity 
for  local  turbulent  characteristics  in  planetary  boundary  layers 
and  the  principle  of  dimensional  analysis  to  perform  a  numerical 
analysis,  and  present  an  estimation  method  for  obtaining  a  vertical 
diffusion  parameter  o2  that  is  applicable  in  the  100  km  range. 

^Comrade  Yen  Pang-liang  has  participated  in  part  of  the  computation 
al  work. 
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II.  DERIVATION  OF  FORMULAS 


According  to  the  hypothesis  proposed  by  Pasquill  [1]  that  the 
local  turbulent  flow  characteristics  are  similar  in  planetary  bound¬ 
ary  layers,  the  average  rate  of  increase  of  i  (the  mean  vertical 
displacement  of  the  diffusing  particle)  is  determined  by  two  local 
parameters  only,  viz.  *•  (standard  deviation  of  vertical  pulsating 

velocity  and  (turbulent  flow  dimension  in  the  vertical  direction) . 
in 

These  two  quantities  are  functions  of  height  and  are,  therefore, 
also  related  to  the  thermodynamic  layer.  One  can  obtain  from  dimen¬ 
sional  analysis 


TT*' 


(1) 


where  f  is  a  function  to  be  determined.  The  rate  of  variation  of 
the  average  horizontal  displacement  x  of  the  corresponding  mass 
point  moving  in  the  same  direction  as  the  wind  is  given  by 

-^-*-*(0  (2) 

Eliminating  dt  from  equations  (1)  and  (2),  one  obtains 

-£-*'(*)  <31 

It  can  be  seen  from  equation  (3)  that  if  #.,#.,  ,  the  distribution 
of  u  (average  wind  speed)  with  height  and  the  functional  form  of  f 
are  known,  then  the  relation  between  z  and  x  can  be  derived  readily 
by  a  numerical  integration. 


1.  Statistical  form  of  the  vertical  vortex  diffusion 
coefficient  K. 


Combining  Taylor's  statistical  theory  and  the  solution  of  Fick's 
diffusion  equation.  We  know  that  <£••*&,  .  After  some  simple  empir¬ 
ical  transformations,  we  can  obtain 

A'-#JL/10  <4> 

tL  is  the  Lagrangian  time  measure.  From  extensive  actual  measure¬ 
ments  [2]  (including  various  different  localities  and  thermodynamic 
layers) ,  we  have 

#•*0-3:'*.  (5) 
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Inserting  equation  (5)  into  equation  (4),  we  obtain 

A'-J^/IS  (6) 

The  above  three  expressions  contain  four  fundamental  quantities: 

e  (power  consumption  rate  of  the  turbulent  flow) ,  K,  o  and  X  . 

w  m 

Given  two  of  these  quantities ,  then  the  other  two  can  be  derived . 

Of  these,  c?w,  Xm  and  c  are  easier  to  measure.  Thus,  one  can  express 
K  in  terms  of  actually  measurable  turbulent  flow  characteristics, 
and  not  have  to  rely  on  hypothetical  mixed  lengths.  Moreover,  the 
above  expression  for  K  is  a  better  representation  of  the  character¬ 
istics  of  the  atmosphere. 

2.  Determination  of  the  function  f. 


As  az  and  K  are  similar  [3],  it  is  easy  to  derive 

d7  _  K 

tlx  «U 

Substituting  equation  (4)  into  equation  (7),  we  obtain 

<h  _  1  ,i- 

tii  io  «  » 

Comparison  of  equations  (8)  and  (3)  shows  that 


It  has  been  given  by  Ito  (1970)  14]  that 

-  * 


(7) 

(8) 

(9) 


Except  for  the  coefficients,  equations  (9)  and  (8)  are  completely 
the  same.  Thus,  it  is  appropriate  to  take  f  to  be  4?-  .  To  impart 
more  generality  to  equation  (8),  we  replace  with  an  undetermined 
parameter  a,  and  regard  a  as  a  function  of  stability.  The  final 
form  of  equation  (8)  becomes 


£.(VU  <10> 

4m  •  '  •  ' 


3.  Selection  of  expressions  for  u,  ow  and  Xn 


Since  the  distribution  with  height  of  ou  and  X.  is  different 

w  TO 

for  different  stability  conditions,  the  expressions  used  below  will 
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be  divided  into  three  groups,  corresponding  to  stable,  neutral  and 
unstable  conditions. 


For  neutral  conditions,  we  adopt  the  relations  obtained  by 
Yohoyama  [5]  in  the  range  of  700-1000  m  above  ground  level: 

•--■■'“•(i— r) 

-t(-4y 

In  the  above  relations,  k  is  von  Karman's  constant,  taken  to  be 
0.35  throughout  the  paper,  h  is  the  thickness  of  the  PBL,  and  u*  is 
a  measure  of  the  velocity  in  the  layer  close  to  the  ground  level. 

He  can  derive  from  equations  (11) ,  (12)  and  (5)  that 


(ID 

(12) 


For  stable  conditions,  we  choose  the  expressions  obtained  by 
Wamser  [2]  from  data  acquired  on  a  300  m  tower: 

*--*•(1— J-)(l.36  +  0.08-£-p)  (14) 

i.-*/(0.36  +  1.3  W*)  ( 1 5 ) 

In  equations  (14)  and  (15),  y  *  jr  is  a  stability  parameter,  and  L 
is  a  measure  of  the  length  in  the  layer  close  to  the  ground  level. 


For  unstable  conditions,  on  the  basis  of  the  results  obtained 
in  references  [11  and  [51- [7],  we  select  the  following  relations 
which  are  applicable  in  the  entire  PBL: 


m-O-tX-'t)* 

i.-* 


(16) 

(17) 


Owing  to  limitations  on  theoretical  understanding  and  available 
experimental  data  of  the  upper  half  of  the  PBL,  there  is  as  yet  not 
a  universally  accepted  expression  for  c?w  and  Xn>  For  instance, 
Pasquill  [1]  assigned  the  number  3.2  to  the  coefficient  in  equation 
(13),  and  pointed  out  that  it  ia  only  applicable  in  the  constant- 
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stress  layer  (usually  100  m  thick) .  Beyond  that,  it  first  increases 
with  height  then  becomes  constant.  From  data  obtained  on  a  tower, 
Wamser  [2]  maintained  that  this  coefficient  is  1.33  at  250  m.  The 
value  2.4  that  we  have  adopted  in  this  paper  was  derived  from  data 
obtained  at  an  even  greater  height  (the  height  attainable  by  verti¬ 
cal  diffusion— 600  m) .  This  value  is  greater  than  1.33  and  smaller 
than  3.2,  and  is  thus  representative  of  an  average  value  for  this 
layer.  Although  the  other  relations  require  further  analysis  and 
verification,  these  were  nevertheless  established  on  the  basis  of 
definite  experimental  data  and  theoretical  analyses.  Using  differ¬ 
ent  expressions  for  different  stability  conditions  is  a  fairly  good 
representation  of  the  problem. 

A  good  description  of  the  wind  speed  in  the  ground  level  layer, 
u(z),  has  already  been  given  in  a  formula  [8).  However,  there  is  as 
yet  not  a  satisfactory  distribution  form  for  the  wind  speed  that  is 
applicable  both  in  the  ground  level  layer  and  in  the  upper  PBL.  In 
recent  studies  on  diffusion  in  the  mesoscale  range,  some  authors  [9] 
have  suggested  the  following  form  for  the  wind  speed  profile: 

«(»)-0(l-«^*  (18> 

In  this  equation,  G  is  the  geostrophic  wind  speed, 

P  is  Coriolis  parameter,  whose  value  is  taken  to 

be  10  s  ,  and  Kq  is  constant  eddy  diffusivity.  In  our  recent  anal¬ 
ysis  of  the  wind  speed  profile  at  a  320  m  tower  [10],  we  found  that 
in  the  320  m  range,  equation  (18)  is  a  better  description  than  the 
power  law  expression.  In  order  to  discuss  the  relation  between  u(z) 
and  the  roughness  and  stability  of  the  substrate,  we  adopted  the 
following  relation  given  by  Tennekes  [11]  for  associating  the  inner 
parameter  uA  and  the  outer  parameters  G,  P  and  z^: 

+  AUriJ*  (19) 

In  the  above  equation,  zQ  is  the  effective  roughness  length  [8],  Q 
is  a  parameter  associated  with  stability,  is  the  Rosby  number 

for  the  ground,  AQ  -  sin  oQ,  and  oQ  is  the  angle  between  the 
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geostrophic  wind  direction  and  the  direction  of  the  wind  in  the 
ground  level. 

Inserting  the  expressions  for  Rg  and  Ag  into  equation  (19) , 
we  obtain,  after  some  operations  and  rearrangement, 

<20> 

In  the  above  equation,  is  the  thickness  of  the  boundary  layer, 

the  value  of  a.Q  varies  between  13  and  35*,  and  to  facilitate  com¬ 
putation  and  at  the  same  time  not  sacrificing  generality,  we  let 
cos a(« i,  .  Thus,  equation  (20)  becomes 

o_^.rlnA.<wj.  (2i) 

III.  EXPRESSION  FOR  FOR  FLAT  UNIFORM  LAND 

Flat  uniform  land  refers  to  the  ideal  condition  where  the 
effective  roughness  length  zQ  =  0.03m.  To  facilitate  comparison 
between  our  results  and  the  available  standard  vertical  diffusion 
types  [13]  (applicable  in  the  10  km  range)  and  to  take  advantage  of 
the  convenience  of  using  the  Gaussian  model,  we  have  adopted  Pas- 
quill's  stability  classification.  We  have  also  converted  the  rela¬ 
tion  between  z  and  x  into  one  between  and  x  with  the  help  of  the 
relation  ».*  az.  The  coefficient  a  is  determined  by  the  vertical 
distribution  of  the  concentration  of  the  pollutant.  Regarding  the 
value  of  a,  Pasquill  in  his  earlier  work  [14]  used  data  from  short 
range  atmospheric  diffusion  experiments  to  obtain  the  average  value 
of  1.3.  Recently,  in  his  study  on  mesoscale  range  vertical  diffu¬ 
sion  [1],  he  pointed  out  that  a  should  be  taken  as  1.25  (correspond¬ 
ing  to  a  normal  distribution  of  concentration  in  the  vertical  direct¬ 
ion)  .  At  the  same  time,  he  pointed  out  that  when  the  index  of  dis¬ 
tribution  y  is  1.5  (2  for  normal  distribution),  a  is  approximately 
1.28,  while  when  the  concentration  has  a  distribution  lying  between 
the  normal  distribution  and  a  power  law  distribution,  a  varies  within 
the  range  1.26-1.42.  This  shows  that,  even  though  the  concentration 
of  the  pollutant  in  the  PBL  does  not  conform  completely  to  a  normal 
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distribution,  our  choice  of  a  =  1.25  for  the  calculations  will  not 
produce  a  large  error. 


I 

i 

I 
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Before  performing  numerical  calculations,  it  is  necessary  to 
make  reasonable  choices  for  the  values  of  the  parameters  encountered 
in  the  calculations.  Our  final  choice  of  these  values  is  given  in 
Table  1. 


The  choice  of  the  values  for  the  thickness  h  of  the  PBL  as  given 
in  the  table  is  based  on  [15]- [17],  height  classification  for  the 
mixing  layer  given  by  Klug  (1969)  [18],  and  the  experimental  results 
obtained  for  the  measurement  of  average  wind  speed  and  temperature 
gradient  on  a  320  m  tower  [19].  L  is  taken  from  the  results  of 
[8]  and  [13] .  The  parameter  1  is  obtained  from  the  relation 

10~*  .  The  different  stability  classes  KQ  are  chosen  from 

the  values  given  by  Draxler  [20]  in  his  study  on  mesoscale  range 
diffusion.  These  values  are  also  given  in  Table  1.  Later,  he  [21] 
also  gave  the  average  values  of  Kg  over  each  of  the  seasons  as  well 
as  the  annual  average.  The  chief  basis  for  the  values  of  Q  is  the 
formula  given  by  Yordanov  [22]: 

eO,)-M+  Mx10->h».is  (22) 

We  have  also  referred  to  [23],  [12]  and  [24]  for  the  selection  of  Q. 
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In  Table  1,  x^  and  z^  are  the  initial  points  for  the  numerical 
integration,  and  have  been  obtained  from  calculations  using  equa¬ 
tions  given  in  [8] .  The  value  of  a  has  been  obtained  by  finding 
by  successive  approximations  that  value  of  a  which,  after  all  the 
other  parameters  have  been  selected,  gives  the  computed  value  of  oz 
that  best  matches  the  curve  given  for  the  range  within  3  km  in  [8], 

The  results  show  that  a  has  smaller  values  for  more  nearly  stable 
atmospheric  conditions.  The  value  for  the  range  of  3-100  km  is 
obtained  by  inserting  the  parameters  given  in  equations  (11) -(21) 
and  Table  1  into  equation  (10)  and  carrying  out  a  numeric  integra¬ 
tion  using  Simpson's  formula. 

After  a  nonlinear  regression  analysis  is  performed  on  the  set 
of  numbers  corresponding  to  and  x  obtained  from  numerical  calcul¬ 
ations,  we  finally  obtain  the  expression  for  *  G(x)  that  is 

applicable  in  the  range  of  0.1-100  km: 

G'(J)-o*7(l  +  «M+"-i>)  (23) 

In  the  above  equation,  a,b,c,d  and  e  are  constants  related  to  stab¬ 
ility.  Table  2  gives  their  values  for  different  stability  conditions. 

IV.  EXPRESSION  FOR  #.  FOR  TERRAINS  WITH  UNIFORM  ROUGHNESS 


In  order  to  analyze  the  expression  for  f#  for  complex  terrains, 
we  adopt  the  concept  of  effective  roughness  length  f8J,  and  use  the 
for  a  terrain  of  uniform  roughness  to  represent  that  for  a  rough 


terrain  with  small  nonuniformity.  Also,  we  write  the  expression 
for  in  the  following  form. 

(24) 

The  G (x)  in  equation  (24)  has  been  given  in  the  previous  section. 

Except  for  Zg,  the  rest  of  the  parameters  remain  the  same.  Repeat¬ 
ing  the  above  computational  process,  and  dividing  the  value  of  •. 
thus  obtained  by  G(x),  we  obtain  values  for  F(x,Zg)  for  different 
stability  conditions,  different  Zg(0.01-5m)  and  different  values  of 
x.  Nonlinear  regression  analysis  then  yields  /309 

•*■(*. *  +C,S  *l  (25) 

The  constants  in  equation  (25)  are  listed  in  Table  3  for  different 
stability  conditions. 
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1 — stability  classification 

It  can  be  seen  from  equation  (25)  that  the  effect  of  Zg  is  pro¬ 
nounced  only  for  short  ranges.  For  very  large  x,  the  right  hand 
side  of  equation  (25)  approaches  1,  and  the  effect  of  local  terrain 
features  of  long  range  roughness  sources  can  be  neglected.  This 
agrees  with  Hanna's  [25]  conclusion.  However,  Hanna  only  used  the 
power  law  form  of  ZgPg  to  express  the  effect  of  Zg  in  *, ,  where  PQ 
varies  between  0.1  and  0.25.  Our  result  is  not  just  a  simple  power 
law  expression,  and  more  closely  represents  the  actual  situation. 

To  clearly  see  the  effect  of  Zg  on  #.  under  different  stability 
conditions,  we  have  listed  in  Table  4  the  values  of  F  for  Zg  =  0.1, 
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Table  4  Ghent*  of  ^ (*,««>  with  (Ubilitiee 
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1  and  5m,  and  x  =  0.1,  1,  10  and  100  km.  From  this,  the  following 
points  become  apparent:  (1)  the  effect  of  roughness  of  terrain 
increases  with  increasing  stability  (A  -*■  E)  .  This  has  been  borne 
out  by  a  large  quantity  of  experimental  data  [26],  and  is  due  to 
the  fact  that  under  unstable  conditions,  oz  is  mainly  controlled 
by  convection  and  not  by  mechanical  turbulent  flow  [25]  .  (2)  if 

zQ  £  1  m,  then  in  the  10-100  km  range  neglecting  the  effect  of  z^ 
will  only  cause  an  error  of  about  25%.  (3)  comparison  of  the  family 

of  curves  for  •,  under  flat  uniform  conditions  and  stability  types 
A-F  with  the  family  of  curves  for  uniform  roughness  shows  that  the 
latter  are  narrower  than  the  former.  In  [27],  a  definite  conclusion 
has  been  drawn  to  this  effect.  It  has  also  been  pointed  out  that 
the  P-G  curves  for  type  A  stability  increases  rapidly  with  x,  and 
are  not  applicable  to  complex  terrains . 

V.  EXAMINATION  OF  RESULTS  AND  COMPARATIVE  ANALYSIS 

In  the  above  two  sections,  we  have  completed  the  calculations 
for  «a  .  The  complete  expressions  are  given  by  equations  (23) -(25) 
and  Tables  2  and  3.  in  this  section,  we  seek  to  answer  the  question 
of  whether  these  expressions  and  choice  of  parameters  are  appropriate 
and  reflect  the  actual  conditions  in  the  atmosphere. 

The  result  given  by  Smith  for  flat  uniform  conditions  is 
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0(i)«a,i»«/(l +«,£*>) 

( 2o 

Equation  (26)  has  been  obtained  by  performing  a  numerical  solution 
on  the  two-dimensional  diffusion  equation  making  use  of  wind  speed 
and  vortex  diffusion  coefficient  profiles  measured  under  different 
stability  conditions. 


The  expression  for  o  used  by  Wendell  in  a  regional  model  is 

z 

0<J)-al*»y(l+«,*+  «»*•)'  <27> 

This  is  a  modification  and  extension  of  Briggs'  [13]  interpolation 
formula.  Although  equations  (23),  (26)  and  (27)  are  not  exactly  the 
same  in  form,  these  agree  in  that  there  is  a  difference  in  diffusi- 
vity  for  short  and  long  ranges,  and  that  the  long  range  diffusivity 
is  smaller.  This  result  is  supported  by  experimental  results  [30] . 
and  shows  that  the  commonly  employed  power  law  form  is  not 

applicable  to  the  mesoscale  range.  For  the  purpose  of  comparison, 
we  have  listed  the  ratios  obtained  by  dividing  equations  (26)  and 
(27)  by  equation  (23)  in  Table  5. 

It  can  be  seen  from  Table  5  that  the  three  models  agree  fairly 
well  in  the  range  of  10-100  km.  For  neutral  conditions,  the  best 
agreement  is  obtained  between  our  results  and  those  of  Smith,  the 
difference  being  within  10%.  For  stable  conditions,  our  results  are 
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smaller  than  Smith's,  with  ratios  lying  below  1.5.  Our  results 
are,  however,  larger  than  those  of  Wendell.  For  unstable  conditions, 
results  obtained  from  our  model  are  slightly  higher  than  those  of 
Smith,  with  ratios  lying  between  0.56  and  1.03.  These  are  smaller 
than  Wendell's.  It  is  remarkable  that  these  three  models  which  have 
been  derived  from  different  angles  of  approach  are  in  such  good 
agreement . 

/311 

We  have  collected  from  related  literature  information  on  oz 
obtained  from  mesoscale  range  measurements  made  at  10  locations  and 
from  estimations  from  the  ground  level  concentration.  There  are  in 
all  72  examples  taken  from  diverse  sources  and  for  different  terrain 
conditions.  Most  of  these  are  for  stable  conditions,  with  48  exam¬ 
ples  for  stability  types  E-F.  The  least  number  of  examples  have  been 
found  for  unstable  conditions  with  only  12  for  stability  types  A-C. 

The  ranges  of  measurements  for  the  latter  are  short  (within  20  km  in 
most  cases) .  An  analysis  of  the  ratio  of  calculated  to  measured 
values  shows  that  72%  of  the  ratios  are  less  than  2,  96%  of  them  are 
less  than  3,  and  only  3%  of  the  data  are  greater  than  3.  The  total 
average  is  0.93  +  0.505.  Comparison  of  Smith's  model  with  the  exper¬ 
imental  data  shows  that  72%  of  the  ratios  are  less  than  2,  81%  of 
them  are  less  than  3  and  19%  of  them  are  greater  than  3.  The  total 
average  in  this  case  is  0.90  +  0.743.  We  conclude  that  there  is 
better  agreement  between  our  model  and  the  experimental  data,  and, 
therefore,  our  model  is  a  better  description  of  the  actual  situation 
than  Smith's  model. 

VI.  CONCLUSION 

The  following  conclusions  can  be  drawn  from  the  numerical  anal¬ 
ysis  given  in  the  above  sections: 

1.  The  estimation  method  of  vertical  diffusion  parameters  in 
the  mesoscale  range  derived  from  the  hypothesis  of  similarity  for 
local  turbulent  characteristics  has  limitations  arising  from  our 
lack  of  complete  understanding  of  the  variation  with  height  of  wind 
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speed  and  the  turbulent  characteristics  in  the  PBL.  Some  of  the 
relations  used  in  this  paper  are  empirical.  Nevertheless,  analysis 
of  the  ratio  of  calculated  to  measured  values  shows  that  there  is 
a  fairly  good  agreement  between  the  two  sets  of  values.  The  average 
ratio  is  0.93  +  0.505.  For  neutral  conditions,  our  results  are 
equivalent  to  those  derived  from  Smith's  model.  For  stable  and 
unstable  conditions,  our  model  is  a  better  description  than  Smith's. 

2.  Comparing  the  variation  of  ».  with  x  in  the  mesoscale  range 

with  that  in  the  10  km  range  we  find  that  this  variation  is  no 

longer  a  simple  power  law  relation,  but  takes  a  complex  functional 

form.  Even  though  the  three  models  that  we  chose  for  comparison 

have  different  forms  of  expression  for  this  variation,  they  all  have 

the  same  tendency,  viz.  the  rate  of  change  of  o_  with  x  becomes 

z 

smaller  with  increasing  distance.  This  result  is  borne  out  by  exper¬ 
imental  data. 

3.  The  effect  of  zQ  increases  with  increasing  stability,  and 
decreases  rapidly  with  Increasing  distance.  For  terrains  with 
roughness  within  z^  ■  1  m,  if  the  effect  of  Zq  is  neglected  in  the 
range  of  10-100  km,  only  a  25%  error  will  result.  However,  the 
importance  of  zQ  becomes  significant  for  ranges  within  10  km. 

In  summary,  although  some  of  the  relations  used  in  this  paper 
are  empirical  in  nature,  And  verification  of  the  results  is  limited 
by  lack  of  adequate  experimental  data,  especially  for  unstable  con¬ 
ditions,  support  for  our  model  is  provided  by  comparison  of  our 
results  with  presently  available  experimental  data  and  satisfactory 
agreement  with  the  results  obtained  from  other  models.  Our  method 
enables  one  to  make  use  of  the  increasing  amount  of  available  data 
on  the  distribution  with  height  of  o .  X  .  e  and  u  to  conveniently 
estimate  the  vertical  diffusion  parameter  without  going  through  the 
complicated  process  of  solving  the  diffusion  equation. 
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